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A  Study  of  Inverse  Methods  for  Processing  of  Radar  Data 


Mohamed  F.  Chouikha,  Ph.D. 


Summary 

This  is  a  technical  report  for  the  project  “A  Study  of  Inverse  Methods  for  Processing  of  Radar  Data” 
supported  by  the  DoD  US  Air  Force  Research  Lab.  This  project  started  on  July  1,  2004  and  continued  until 
August  31,  2005.  The  goal  of  this  project  is  to  investigate  the  possibility  of  new  inversion  algorithms  for 
radar  image  processing  to  improve  signal  quality  and  reduce  the  effects  of  clutter  based  on  study  of  known 
geophysical  inversion  algorithms. 

1.  Introduction 

The  scattering  of  energy  can  be  considered  on  a  point  by  point  basis,  each  diffraction  point  contributing  a 
part  of  the  total  signal.  Kirchhoff  methods  pioneered  by  Bleistein  and  Cohen  at  the  Colorado  School  of 
Mines  [6]  used  an  asymptotic  integral  method  to  construct  a  migration  image  of  the  seismic  data.  By 
reduction  for  the  number  of  integrations  they  were  able  to  make  this  method  computationally  feasible.  The 
method  needs  accurate  travel  times  for  the  primary  wave  front  to  every  subsurface  point  from  each  source 
receiver  location.  Both  ray  tracing  and  eikonal  schemes  have  been  used  to  computer  these  travel  times.  A 
by-product  of  their  scheme  is  a  local  estimate  of  the  parameters  beneath  the  reflector  surfaces  constructed 
from  two  different  integral  approximations. 

Claerbout  [9]  and  his  Stanford  Exploration  Project  took  the  approach  of  downward  continuing  the  surface 
data  using  a  paraxial  approximation  and  its  extensions  to  larger  angles.  Their  further  work  in  this  area  has 
developed  angle  versus  aperture  imaging  algorithms  using  the  one  way  wave  equation. 

The  use  of  angle  dependent  plane  waves  for  seismic  imaging  was  invented  at  the  Amoco  Research  Center 
by  Dan  Whitmore.  He  was  able  to  show  in  his  PhD  thesis  that  using  a  few  plane  wave  angles  and 
converting  the  time/distance  data  into  ray  parameter  and  time  or  tau-p  data  that  very  accurate  images  can  be 
developed.  Gazdag  (1978)  [7]  and  then  Stoffa  et  al.  (1990)  [4]  made  use  of  phase  shift  and  interpolation 
phase  shift  schemes  to  develop  very  good  migration  algorithms,  again  using  the  one  way  wave  equation. 

All  of  these  methods  make  transformations  to  the  wave  equation.  Another  computer  expensive  but  very 
successful  method  uses  the  approach  of  running  the  wave  equation  backward  in  time.  This  uses  the  inverse 
of  the  modeling  equation  for  imaging.  First  published  by  Hemon  [1 1],  it  was  popularized  by  Whitmore  [13] 
after  Kelly  and  team  [12]  demonstrated  successful  acoustic  modeling  codes.  Currently  for  real  three- 
dimensional  elastic  imaging,  prestack  reverse  time  methods  [  1 0]  are  the  most  accurate  and  are  competitive 
in  cost. 

A  distinct  characteristic  of  geological  interface  surfaces  is  the  relatively  long  specular  surfaces  in 
relationship  to  the  surface  source/receiver  geometry.  To  extract  subsurface  features  like  river  channels  and 
fractures  from  such  regular  structures  has  posed  a  significant  challenge.  Several  researchers  in  the  seismic 
arena  have  constructed  coherency  measures  to  highlight  these  relatively  short  features.  The  eigenvalue 
method  of  Gerztenkom  is  typical  of  these  seismic  coherency  methods.  Within  the  Fresnel  zone  the 
resolution  of  good  quality  seismic  data  is  thought  to  be  on  the  order  of  a  quarter  of  a  wavelength.  With  the 
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much  higher  frequencies  of  radar  and  at  the  resolution  needed  for  imaging  vehicles  the  expected  number  of 
wavelengths  in  a  typical  data  volume  will  be  much  larger  than  in  seismic.  This  will  impose  further 
restrictions  on  the  use  of  these  seismic  based  methods,  but  modeling  studies  will  demonstrate  with  these 
algorithms  their  relative  merits. 

Signal  point  diffractions  have  known  responses  and  by  using  time  migration  or  depth  migration  these  can 
be  collapsed  to  the  scatter  point.  The  use  of  numerous  points  in  modeling  along  the  boundaries  to  act  as 
absorber  or  sponges  is  the  basis  for  a  numerical  boundary  condition.  This  field  of  scatter  points  is  regular 
and  recent  work  shows  clearly  the  effects  of  dense  scattering  points  [14].  How  this  applies  to  radar  clutter 
scattering  is  not  yet  known  but  the  seismic  response  is  dramatic  for  both  acoustic  and  elastic  wave 
equations.  Further,  we  have  developed  a  mathematical  basis  for  the  numerical  method,  creating  a  link  to 
the  previous  analytic  methods.  Essentially,  the  sponge  boundary  condition  has  the  features  of  the  angle 
dependent  analytic  methods  with  a  systematic  scheme  for  specifying  the  angles. 

In  this  report,  we  focus  on  reviewing  a  number  of  traditional  seismic  algorithms,  especially  the  finite- 
difference  reverse-time  migration  algorithm  and  the  phase  shift  migration  algorithm.  We  also  perform 
preliminary  experiments  on  synthetic  seismic  data  using  these  algorithms  and  analyze  their  implications  in 
radar  signals  processing.  We  further  discuss  their  merits  and  shortcomings,  and  provide  a  plan  for  future 
research  of  this  project. 

2.  Algorithms  review 

2,1  Finite-difference  reverse  time  migration  algorithm 

This  method  is  one  of  the  depth  migration  methods  that  use  finite  difference  solutions  to  the  acoustic  or 
elastic  wave  equation  with  the  time  reversal  of  the  seismic  traces  being  employed  as  time  varying  sources  at 
the  surface  boundary.  Migration  is  a  technique  to  position  seismic  reflections  in  their  correct  subsurface 
location  based  on  recorded  seismic  traces  on  the  surface.  A  reflector  can  be  treated  as  a  large  set  of  point- 
diffractors  in  the  subsurface  when  seismic  waves  are  refracted.  This  time  reversal  allows  the  wave  field  to  be 
back  propagated  to  the  depth  point  where  the  reflections  originated.  Finite-difference  reverse  time  migration 
algorithm  is  the  most  general  method  of  all  depth  migration  codes  although  it  is  also  the  most  computer 
intensive.  Using  the  complete  wave  equation  allows  energy  to  move  in  all  directions  during  the  back 
propagation  process.  The  one  way  wave  equation  methods  only  allow  propagation  into  the  model  from  the 
boundary  surface.  If  energy  is  “multiply  reflected”,  bouncing  back  and  forth  between  impedance  surface 
changes  the  complete  wave  formulation  allows  for  these  effects. 

The  basis  for  this  migration  of  the  stacked  seismic  time  section  is  the  exploding  reflector  model  [  1  ].  For 
zero-offset  case  (coincident  source-receiver  situation),  the  response  for  a  wave  field  propagating  up  from 
the  reflector  at  half-velocity  of  the  medium  is  considered  equivalent  to  the  a  wave  fields  propagating 
downwards  from  the  source  to  the  reflector  and  then  back  to  the  receiver  at  the  medium  velocity.  So  the 
exploding  reflector  model  assumes  that  if  the  downward  wave  travels  the  same  path  as  the  reflected  upward 
wave,  then  only  one  wave  propagation  is  needed.  The  upward  wave  will  be  chosen  for  modeling  because  it 
is  the  one  that  is  recorded  at  the  surface,  and  the  migration  will  attempt  to  reconstruct  the  depth  image  of 
the  exploding  reflector  wave  field  from  the  knowledge  of  the  surface  recorded  wave  fields. 

A  migration  problem  can  then  be  viewed  as  an  attempt  to  use  data  recorded  at  earth  surface  z=0  in  order  to 
reconstruct  the  depth  image  of  the  exploding  reflector  at  t=0.  In  other  words,  migration  attempts  to 
reconstruct  the  exploding  reflector  wave  field  u(x,y,z,t=0)  from  the  knowledge  of  earth  surface  recorded 
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wave  field  u(x,y,z=0,t)  [3].  Good  results  with  reverse  time  migration  depend  on  general  algorithm,  good 
input  data,  and  an  accurate  velocity  model. 

(a)  Forward  modeling  (for  data  generation) 

The  zero  offset  data  satisfy  the  wave  equation 


dt2 


c2(x,z) 


d2V 

dx2 


(1) 


where  c(x,  z)  is  the  velocity  function  (acoustic  sound  speed) 

^f(x,z,t)  is  the  pressure  wave  fields 

(x,z,t)  are  surface,  depth,  and  time  coordinates. 

This  equation  provides  a  relationship  between  the  second  temporal  derivatives  and  the  second  spatial 
derivatives.  The  finite  difference  approximation  of  the  wave  equation  (1)  is 

_  iyvN  4.  VUN~l 

-  (2) 

Solving  for  the  new  time  N+l,  we  get 


'PA'+1  =2vP"-'f lN~'+St2c2 


d2vr 

dx2 


+  src(t) 


(3) 


where  src  (t)  is  a  source  function  which  helps  in  generating  artificial  seismic  pressure  waves  and  in 
recording  their  reflections  on  the  surface.  The  result  from  this  step  will  be  the  synthetic  data  (time  section). 
An  example  source  function  is  given  below  which  is  plotted  in  Figure  1 . 


src(t)  =  sin(2^/nuxr) 


e(-a,2) 


P 
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(b)  Reverse  Time  Imaging 

The  seismic  data  are  introduced  to  the  wave  equation  as  boundary  condition  (data  would  propagate 
into  earth  from  the  subsurface  receiver  position),  and  the  wave  equation  is  then  formulated  as  a  backward  in 
time  propagation  until  time  reaches  zero  (imaging  condition). 

Using  the  exact  same  wave  equation  but  now  solving  for  'F  v  /,  we  get 


ipw-i  =2'¥N  -V"*'  +St2c2 


d2xV 
dx 2 


+  - 


a2T 

dz2 


+  B.Ck 


J  N 


(4) 


2.2  Phase  shift  migration  algorithm 

The  phase  shift  method  of  migration  is  based  on  the  numerical  solution  of  the  wave  equation  in  the 
frequency  domain  with  the  initial  conditions  defined  by  the  zero  offset  seismic  section.  The  zero  offset  data 
are  the  wave  fields  {ft^x,  z=0,  t)}  measured  at  some  specific  depth  from  the  surface  of  the  earth.  This 
method  begins  with  the  Fourier  transform  of  the  data  set.  This  transformed  data  is  then  extrapolated 
downwards  and  subsequently  evaluated  at  t=0.  This  method  is  a  one  way  transform  of  the  wave  equation. 

In  case  of  horizontal  layered  velocity,  migration  is  defined  by  a  set  of  independent  ordinary  differential 
equations  in  the  k-w  domain,  and  the  wave  components  are  extrapolated  downwards  by  rotating  their 
phases  [7].  In  case  of  lateral  velocity  variations,  the  wave  field  is  extrapolated  by  phase  shift  methods  using 
n  laterally  uniform  velocities.  The  intermediate  result  is  n  reference  wave  fields.  Then  the  actual  wave  field 
is  computed  via  interpolation  from  the  reference  wave  fields  [8]. 

(a)  Horizontal  Velocity  Variation 

The  pressure  wave  fields  Vfjc,  z=0,  t)  (zero  offset  data)  recorded  at  the  surface  satisfy  the  wave 
equation 

d2xi>  _  1 

dz2  ~  v2 

The  process  of  migration  is  reconstructing,  at  the  reflectors,  the  signals  that  produce  the  seismic  section 
when  propagated  to  the  surface,  in  other  words  trying  to  evaluate  ¥fjc,  z,  t=0 ).  These  values  are  obtained  by 
solving  equation  (5)  with  initial  values  ffx,  z-0,  t). 

To  solve  for  ffx,  z,  t=0),  let  ffx,  z,  t)  be  represented  by  Fourier  series 

<F(r,z,/)  =  ^^'f(^,z,  go)  exp  [i  ( kxx  +  cot )]  (6) 

kx  (0 

where  kx  is  the  midpoint  wave  number,  and  co  is  the  frequency. 

Substituting  (6)  to  (5)  gives 


a2vF 

dt 2 


d2vF 

dx2 


(5) 
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(7) 


d2vr 

dz 2 


=  -k'¥ 


The  analytical  solution  to  equation  (7)  is 

X\’(kx,z  +  Az,a>)  = 'F(&t,z,<y)exp[/A:,Az]  (8) 

w 

where  k,  =  ±— 
v 


'-f-T 

v  w ; 


72 


It  should  be  noted  that  the  two  values  of  k,  correspond  to  forward  time  (negative  value)  and  reverse  time 
(positive  value)  wave  propagation.  Downward  extrapolation  requires  that  A z  be  greater  than  zero.  Since 
downward  extrapolation  of  the  recorded  seismic  data  is  an  inverse  problem,  only  the  positive  value  of  k,  is 
used  in  the  solution. 


The  final  expression  for  the  wave  extrapolation  is  thus  given  by 


w 


Xi'(kx,z  +  Az,  w)  =  xV{k(,z,  w)exp<z — 

v 


w  ) 


iK 


Az 


(9) 


This  is  the  solution  to  the  wave  equation 

dV 


'to' 


dz 


=  i 


\vj 


1- 


feY 

<  a>  7 


1 X 


¥ 


(10) 


Equation  (9)  is  considered  the  basis  of  the  phase  shift  method,  and  equation  (10)  is  the  exact  extrapolation 
equation  for  the  case  of  depth  velocity  variation  with  v=v  (z). 

Finally,  substituting  equation  (9)  into  the  Fourier  series  (6)  and  evaluating  at  t=0  will  provide  the  solution 
to  the  migration  problem  for  the  depth  variable  velocity  case. 

'F(x,z  + A z,t  =  0)  =  ^^4/(^,z,<w)exp[/(A:xx  +  A:TAz)]  (1 1) 

kx  to 

Equation  (6)  can  also  be  written  as 

^(x,  z,  t  =  0)  =  ^  £  ^(^ ,  z  =  0,  ai)  exp  [/  (£,*  +  k2z )]  (1 2) 

kx  o) 

where  vF(x,z  =  0,t),  the  coefficient  of  the  series,  represents  the  Fourier  transform  of  the  seismic  section. 
(b)  Lateral  Velocity  Variations 

In  this  case  the  solution  expressed  by  equation  (9)  is  no  longer  valid  for  the  fields  with  lateral 
velocity  variations,  and  the  square  root  of  equation  (10)  is  expanded  into  truncated  series  (i.e.  Taylor 
series). 
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3.  Preliminary  experiments 

We  have  performed  some  preliminary  experiments  using  the  code  developed  for  reverse  time 
migration  algorithm  and  the  phase  shift  migration  algorithm  as  provided  in  the  open  source  Seismic  Unix 
package. 

3.1  Reverse  time  migration  algorithm 

This  algorithm  is  illustrated  below  with  synthetic  data.  A  layered  salt  synthetic  model  is  shown  in 
Figure  2  and  the  corresponding  migrated  image  using  reverse  time  migration  algorithm  is  shown  in  Figure 
3. 


Depth  (z) 


Surface  coordinate  x 


Figure  2  A  salt  model  for  synthetic  data  generation 


Figure  3  Migrated  image  for  model  in  Figure  2 

Figure  4  shows  the  input  velocity  model  which  consists  of  three  reflector  segments,  a  horizontal 
reflector  and  the  other  two  reflectors  with  dips  45  and  60  degrees.  The  velocity  model  is  a  simple  constant 
model  with  v=1500  m/sec,  and  the  model  dimension  is  400x200  with  grid  spacing  Ax  =25,  Az=25. 
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Surface  coordinate  x 


Figure  4  Input  velocity  model 


The  time  section  or  the  data  resulting  from  the  forward  modeling,  and  the  migrated  image  for  the  model  are 
shown  in  Figures  5  and  6,  respectively. 


Surface  coordinate  (x) 


Time  (t) 


Figure  5  Zero-offset  seismic  section  on  the  input  model 
Surface  coordinate  (x) 
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As  can  be  seen  from  the  above  figures  it  is  noticed  that  the  dips  and  the  location  of  the  reflectors  and  layers 
agree  with  those  in  the  original  models. 

4.  Discussion 

The  scattering  of  waves  by  clutter  causes  noise-like  effects  mixed  with  coherent  signal.  Because  of  the 
multiple  like  processes  which  are  familiar  in  geophysical  processing  of  seismic  data  we  believe  that  using 
the  complete  wave  equation  should  be  our  first  attempt  at  modeling  and  migration. 

An  exploding  reflector  and  reverse  time  migration  (RTM)  code  was  developed  for  the  2D  acoustic  wave 
equation.  This  code  allows  for  a  heterogeneous  velocity  model,  a  general  gridded  function  of  X  and  Z. 
Geophysics  uses  X  for  the  horizontal  axis  and  Z  for  the  depth,  positive  downward.  Our  first  experiment  is 
to  use  a  blocky  earth  model  to  generate  a  zero  offset  (X,T)  data  set.  Then  using  the  reverse  time  migration 
code  as  described  above  we  generated  the  migrated  image  of  the  original  model.  Of  course  if  the 
parameters  are  known  exactly  this  should  be  an  excellent  recovery  of  the  model  as  we  demonstrate.  The 
next  experiment  is  to  use  a  constant  velocity  background  and  a  set  of  finite  reflectors  placed  at  different 
angles.  These  are  then  used  as  the  source  locations  for  the  exploding  reflector  modeling  phase.  The  RTM 
of  this  data  gives  an  excellent  result. 

We  plan  on  creating  point  diffractors  in  the  region  around  these  finite  reflectors  to  simulate  the  effect  of 
clutter.  Currently,  we  are  investigating  small  variations  of  the  velocity  model  around  a  central  mean  value. 
These  velocity  gradients  are  intended  to  simulate  the  clutter  velocity  shift,  the  reduction  in  velocity  due  to 
scattering.  These  gradients  are  seen  in  ground  penetrating  radar  (GPR)  data  and  significantly  better  images 
are  obtained  in  the  GPR  images  when  small  changes  in  velocity  are  considered. 

The  phase  shift  methods  for  migration  are  suitable  for  velocity  models  which  vary  in  only  one  dimension. 
These  v(z)  models  are  simple  to  implement  in  the  method  as  each  step  downward  requires  a  z  increment. 
These  methods  are  being  tested  with  the  same  exploding  reflector  zero-offset  data  that  is  being  used  for  the 
RTM  process. 

5.  Conclusions 

This  is  a  technical  report  for  the  project  “A  Study  of  Inverse  Methods  for  Processing  of  Radar  Data” 
supported  by  the  Air  Force  Research  Lab.  This  project  started  on  July  1,  2004  and  completed  August  31, 
2005.  The  goal  of  this  project  is  to  investigate  the  possibility  of  new  inversion  algorithms  for  radar  image 
processing  to  improve  signal  quality  and  reduce  the  effects  of  clutter  based  on  study  of  known  geophysical 
inversion  algorithms.  This  technical  report  summarized  what  has  been  done  concerning  algorithm  review 
and  testing  during  this  period  of  time. 
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